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Maximally localized Wannier functions are widely used in electronic structure theory for analyses 
of bonding, electric polarization, orbital magnetization, and for interpolation. The state of the 
art method for their construction is based on the method of Marzari and Vanderbilt. One of the 
practical difficulties of this method is guessing functions (initial projections) that approximate the 
final Wannier functions. Here we present an approach based on optimized projection functions that 
can construct maximally localized Wannier functions without a guess. We describe and demonstrate 
this approach on several realistic examples. 


I. INTRODUCTION AND MOTIVATION 

Within the quasiparticle approximation, the electronic 
states of a crystal can be described in terms of single¬ 
particle Bloch functions if) n k(r). These functions are 
eigenstates of the crystal Hamiltonian, and can be labeled 
by their band index n and crystal momentum k. Wannier 
functions (WFs) provide an alternative representation in 
which an entire band of electrons is described by a single 
function |Rn) localized in or near the unit cell labeled 
by the lattice vector R. In their simplest form, WFs are 
obtained from the Bloch functions via the Fourier trans¬ 
formation 

l Rn ) = TWA / dke“* k ' R \ip n k) , (1) 

(2tt) J bz 


the WFs |Rn) in r space will depend on the smoothness 
of the gauge in k space. If the V’nk(r) are chosen 
with random overall complex phases [which often hap¬ 
pens if ip n k(r) are acquired numerically by diagonalizing 
a k-dependent Hamiltonian matrix separately for each k 
point] the WFs obtained from Eq. (1) need not be local¬ 
ized. However, if matrices are chosen so that Bloch 
states are smooth in k space (smooth gauge), the corre¬ 
sponding WFs will be localized in r space. 

The idea of maximally localized Wannier functions and 
a procedure for obtaining them from a set of composite 
bands was introduced by Marzari and Vanderbilt' for 
isolated bands and later extended to the case of entan¬ 
gled bands. 8 Maximally localized Wannier functions are 
constructed by choosing a gauge Umn for Eq. (2) that 
minimizes the spread functional 


where V is the volume of the real-space primitive cell. 
The definition of WFs is not unique because there is a 
gauge freedom in the right-hand side of Eq. (1). Namely, 
at each k point and for each n, one can change the overall 
phase of the Bloch state I^Vsk}- In fact, one often consid¬ 
ers an even more general gauge choice which allows an 
arbitrary unitary transformation of a set of N bands at 
each k point, 


IV’nk) ’ ( 2 ) 

m 

We focus here on the case when these N bands are iso¬ 
lated from the rest. The choice of gauge is now expressed 
through a k-dependent N xN unitary matrix . 

When Wannier functions are localized in real space 
they have a wide use in the electronic structure commu¬ 
nity. An extensive review of maximally localized Wannier 
functions (MLWFs) and their properties and applications 
can be found in Ref. 1. For example, they have been used 
in the description of electronic polarization 2 and orbital 
magnetization, in addition to being used for interpolation 
of band structures and matrix elements 3-5 and electron 
transport calculations. 6 

For this reason, one often uses the gauge freedom 
so that the corresponding WFs are localized. As a general 
consequence of the Fourier transform, the localization of 


n = £[«„-3]. (3) 

n 

where 

( r2 ) n = ( 0n \ r2 \ 0n ) ’ ( 4 ) 

r n = (0n|r|0n) . (5) 


Here the spread functional fi is written in terms of the 
Wannier functions |0?r). Usually there exists a global 
minimum of Ul corresponding to a unique choice of 
(up to translation of the WFs and their overall complex 
phase), but in some cases there are multiple solutions.' 

Using the general form of Eq. (1) including the u *- k - ) 
matrix in Eq. (2), the spread can be recast in terms of 
the Bloch states. More specifically, Ul can be expressed 
only as a function of the overlaps of the periodic parts of 
the Bloch functions at neighboring k points k and k + b, 

m if h) = («ik|«jk+b) • (6) 

See Appendix A for an explicit definition of S2 in terms 
of 


n = n 



( 7 ) 


Here we only note that spread H can be decomposed into 
three parts: the invariant part, which does not depend 
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on the gauge, and the diagonal part and the off-diagonal 
parts which do, 


ci 


m 


(k,b) 
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The procedure for minimizing Cl, outlined in Refs. 7 
and 8, is implemented in the Wannier90 code 9 and 
has become the standard method for obtaining localized 
WFs. A notable drawback in the standard approach that 
we address in this manuscript is that one often needs to 
provide a good initial guess of the MLWFs to find the 
global minimum of Cl. In this work we demonstrate a 
modified procedure, in which localized Wannier functions 
are constructed as a linear combination of physically 
based atom-centered orbitals without requiring an ini¬ 
tial guess, as in the standard approach. 1 This is achieved 
by finding optimal projection functions (OPFs) so that 
the resulting Wannier functions obtained via projection 
(as in Sec. II A) are as localized as possible. This OPF 
method could, for example, be used in constructing mate¬ 
rial properties databases, such as the database of the Ma¬ 
terials Project, 10 by providing a simple localized Hamil¬ 
tonian that could serve as a descriptor for the electronic 
structure of a material. We present the theoretical ap¬ 
proach and numerical methods in Sec. II and III. Several 
realistic materials are investigated in Secs. IV to illus¬ 
trate our approach for constructing localized WFs. 

Schemes beyond the standard implementation '~ 9 have 
been developed by others to improve the construction 
of MLWFs and their properties. The inclusion of un¬ 
occupied anti-bonding states has been shown 11 to give 
more localized Wannier functions, but at the expense of 
a chemical picture of the occupied states. Additionally, 
constraints on the matrices can be imposed in order 
to construct localized Wannier functions that possess all 
the space group symmetries of the crystal. 12 


II. STANDARD APPROACH 


Here we summarize the main result of Ref. 7 for a two- 
step construction of maximally localized Wannier func¬ 
tions. In the first step of minimizing the spread func¬ 
tional one needs to guess orbitals g :] (r) with roughly 
the same orbital characters and real-space location fj as 
the target MLWFs. This choice is often done based on 
an intuitive understanding of the band structure of the 
crystal under investigation. Given a choice of gj( r) close 
to target WFs, one constructs the gauge for which spread 
functional Cl is near its global minimum (better choices 
give Cl closer to the global minimum). In the second step, 
this initial gauge choice is iteratively optimized until Cl 
reaches a global minimum. In practice, the second step 
usually reduces the spread fl only by 20% or less. 


A. First step 

Now we describe the first step of this procedure in the 
simple case of a single band of states ipk ( r )- Given a 
localized function g(r) approximating the target MLWF 
at the origin, we first project it onto the Bloch state '0k(r) 
at each k 


fl (k) = (V>k| g) • (9) 

Now we rotate the phase of Bloch state V'k(r) so that the 
relative phase of rotated Bloch state and g( r) is zero for 
all k points, 


|W^a (k) (a (k) *a (k) ) V V). (10) 

It is easy to check that if the initial guess g( r) were a true 
target MLWF, inserting these rotated Bloch states into 
Eq. (1) would give back the target MLWF. However, since 
g(r) is only an approximation, the spread Cl of the rotated 
Bloch states is not exactly at the global minimum. For a 
good guess g{ r), however, the spread should be close to 
the global minimum. 

Following the procedure in the case of a single band, 
we now generalize it to the case of N composite bands. 
First, we choose a set of N localized orbitals gj( r) that 
are approximately equal to the N target MLWFs, 


1 9j) ~ |0 j) ■ (11) 


Here we choose for convenience \gj) to be close to the 
MLWFs near the origin (R = 0), but in principle any 
other R can be chosen. 

Next we compute the overlap between all N Bloch 
bands and N initial guesses for the WFs, 

a- k) = (Vwlffj) • (12) 


Unlike the case of a single band, a (k) for an isolated group 
of N Bloch bands is a N x N matrix, so that Eq. (10) 
generalizes to 


iV’ik) u p ] l^ik) = a i l ( at °) 

j ji 


- 1/2 


IV’jk) • 


(13) 

Here, to simplify notation, we suppress the k dependence 
of a. The inverse square root on the right-hand side is the 
matrix square root of (o'a) . For further simplification, 
define u x , for an arbitrary matrix x, as 


u x = x(x^x) . 


(14) 


The matrix u x is unitary by construction. In fact, it is 
the closest unitary approximation of x. 

In practice, the unitary matrix u x is obtained via the 
singular value decomposition (SVD) of x. If x = zdv 
is a SVD with z and v unitary, and d diagonal, then 
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u x is simply zv. Given this notation, the earlier gauge 
transformation from Eq. (13) now reads 


iV’ik) -> 

3 



(15) 


As in the case of a single band, if trial orbitals gj( r) were 

(k) 

chosen close enough to the target MLWF, the gauge u\ 
from Eq. (15) will by definition give a spread close to the 
global minimum, 


n = n 


u (k)t TO ( k , b ) u (k+ b ) 


(16) 


Here we implicitly wrote in terms of overlap matrices 
m( k,b ) and used gauge transformation of Bloch states 
from Eq. (2) to get transformation of the overlap matrix 
m (k,b) defined in Eq. (6). 


B. Second step 


The initial gauge Ua J can be further improved in the 
second step by rotating, at each k point, the gauge 
from ui k) to u^v l ' k ' > with an appropriate choice of k- 
dependent matrices The spread functional Q is 

minimized using the method of steepest descent. The 
gradient is determined by calculating the derivative of 
the spread with respect to the unitary matrices ?/ k ) and 
then following the path along the direction which mini¬ 
mizes fi. Written more formally, the second step of the 
standard procedure finds a set of unitary NxN matrices 


€ U(N,N), one for each k, 


that 


minimize 


n 


u (k)t 4 k)t m (k ’ b) ui k+b) u (k+b) 


(17) 


(18) 


Quite generally, the global minimization of a function 
using the steepest descent algorithm is bound to work 
well when one starts near the global minimum. Otherwise 
it is quite possible for the algorithm to get stuck in a 
local minimum. In other words, the second step of the 
procedure will arrive at the true MLWFs as long as the 
initial guesses gj (r) in the first step are close enough. It is 
this issue that we aim to address in this manuscript: how 
to automatically construct a gauge that is guaranteed to 
be close to the global minimum. 


III. ALTERNATIVE APPROACH 


combination of hj. In other words, the space spanned by 
hj must approximately contain, as a subset, the space 
spanned by the MLWFs near the origin, 

Span(|/ij)) D Span(|On)). (19) 

The requirement on | hj) is significantly less restrictive 
than that on \gj) in the standard approach. In fact, the 
requirement Eq. (19) should be rather easily satisfied. 
Since we expect MLWFs to be linear combinations of 
atomiclike valence electrons, we can simply choose hj 
to be a set of atom-centered atomic orbitals for each 
atom in the crystal basis and for each relevant atomi¬ 
clike orbital in the valence (some combination of s, p , d, 
/ atomic orbitals, depending on the valence). If nomi¬ 
nal valence atomiclike orbitals are not enough to satisfy 
Eq. (19) (which might happen for example in material 
under extreme pressure), one can always include atomic 
orbitals with higher radial and orbital quantum numbers. 

In the case of covalently bonded materials, a specific 
target MLWF might have its center on a covalent bond at 
the edges of the primitive unit cell. If this is the case, then 
we can expand the set hj(r) by including the periodic 
images of a few atoms in the crystal basis, so that in the 
end, for each unique covalent bond, both atoms forming 
the bond are included in hj( r). 

Since the functions hj satisfy Eq. (19), it is possible to 
approximate the MLWFs as linear combinations of hj. 
Formally, it is possible to find a semiunitary rectangular 
MxN matrix W such that the functions 

M 

Ifc-HEWylfc) (2°) 

i =1 

are close to the target MLWFs. (Since W is rectangular, 
the N functions gj are linear combinations of M functions 
hj.) Thus obtaining approximate MLWFs is equivalent 
to finding the matrix W. We shall call these gj optimized 
projection functions (OPFs). 

To measure the closeness of gj to the target MLWFs, 
we need to express spread in terms of W. Therefore, 
we first need a projection of gj into Bloch states. Since 
g depends on W, it is more convenient to first project 
hj onto the Bloch states, yielding the NxM projection 
matrix 

4 k) = <^*> • ( 21 ) 

Given A (k ) we can compute the overlap matrix between 
the gj and the Bloch states, 

M 

4^ = (V’ikl gj) = {^ik\hl) Wij , (22) 

i=i 


In our approach, instead of choosing N functions g 3 (r) 
that are close to the N target MLWFs, we start with a 
larger set of M functions (M >N) labeled hj( r). These 
functions hj will be chosen so that any MLWF near the 
origin (R = 0) can approximately be written as a linear 


or, in short, 

a ( k)=A (k )w. (23) 

Here we adopted the convention that small (NxN) square 
matrices are written with lower-case Latin letters, while 
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rectangular (TV x M or M x TV) or large square matrices 
(M xM) are denoted by upper-case Latin letters. 

Now we are ready to express fi in terms of W. Com¬ 
bining Eq. (16) and Eq. (23) yields 


n = n 


(k)f (k.b) (k+b) 


l AW 


l AW 


(24) 


To draw comparison with Eqs. (17) and (18), in our ap¬ 
proach the process of constructing MLWFs is equivalent 
to finding 


W € U(M,N), a single matrix 


that 


minimizes 


n 


(k 2t TO ( k ’ b )J k _± b ) 


w 


‘Aff 


(25) 


(26) 


Once the W which minimizes Eq. (26) is found, we use 
the matrices to rotate Bloch states at each k point 
into a smooth gauge. In most of the concrete cases stud¬ 
ied, the spread of the Wannier functions corresponding 
to this gauge is within 1% of the global minimum (this 
is discussed further in Sec. IV) and therefore there is no 
need to improve the gauge further. However, in principle 
one could run the second step of the standard procedure 
to bring spread to its true global minimum and thus ob¬ 
tain maximally localized Wannier functions. 

Now we will compare our approach to the standard 
method in more detail, outlining both the advantages 
and disadvantages of our approach. We also discuss the 
approximations that are made to implement an algorithm 
to construct the W matrix. 


A. Comparison to the standard approach 

The procedure for constructing MLWFs by generating 
OPFs [Eqs. (25) and (26)] has several advantages com¬ 
pared to the standard procedure [Eqs. (17) and (18)]. 
First, OPF construction is given by a single matrix W, 
instead of a set of ib k ) matrices, one at each k point. 
For this reason, as will be shown in Sec. IIIB, one can 
more directly solve Eq. (26) without using the method 
of steepest descent; rather, an iterative procedure is used 
to construct IE as a product of large unitary transfor¬ 
mations (Givens rotations). Therefore, this procedure is 
less likely to get stuck in a local minimum. The second 
advantage of OPF construction is that the W matrix it¬ 
self has a lot of chemical information encoded in it. For 
example, one can see directly from W the contribution 
of the various atomic orbitals to each OPF and thus the 
corresponding Wannier functions. We discuss this point 
on concrete examples in Sec. IV. Third, the use of a single 
matrix might make it easier to impose constraints such 
as crystal symmetry. 

There are however some disadvantages to the OPF con¬ 
struction approach. First, the spread II in Eq. (26) de¬ 
pends nonlinearly on W since it appears under the matrix 


(k) 

inverse square root in u aw- I n f &c tj Taylor expansion of 
the inverse square root leads to a power series in all pos¬ 
itive integer powers of W. Second, since we do not want 
to rely on a steepest decent method, minimization of the 
diagonal part of the spread [IId in Eq. (8)] becomes non¬ 
trivial. 

In the following section we introduce two simplifica¬ 
tions to Eq. (26) which deal with these two disadvantages 
of OPF and allow for an efficient numerical construction 
of OPFs in all the cases studied. 


B. Simplifications 

The following two subsections describe two simplifica¬ 
tions that turn minimization of Eq. (26) into a numeri¬ 
cally efficient form. 


1. Linearizing ilaw 

The first simplification in minimizing the spread IT 
from Eq. (26) is to expand it to the leading order in W. 
Explicitly writing uaw in terms of its definition [Eq. (14)] 
and ignoring k index for the moment, 

uaw = AW (W t A*AW)~ 112 . (27) 

For W which minimizes Eq. (26) we expect 

IVU (k)t A (k) fE « I N , (28) 

for all k since the OPFs approximately overspan the 
space of MLWFs (In is the NxN identity matrix). There¬ 
fore, at least near the optimal value of W, we are justified 
in Taylor expanding ilaw around W^A^AW close to the 
identity (In), 


uaw = AW 


I N -^(wUUW-I N ) + ... 


(29) 


Therefore, to lowest order, uaw ~ AW. Restoring uni- 
tarity we can replace A with Ua , thus obtaining a unitary 
approximation to Uaw > 


uaw » U A W. (30) 


Here Ua has been constructed according to the Lowdin 
orthonormalization procedure given by Eq. (14). We fol¬ 
low here the notation we introduced earlier so that Ua 
with upper case U is a rectangular NxM matrix (while 
uaw with lowercase a is a square TV x TV matrix). We 
also note here that Eq. (30) is exact if W were a square 
matrix. 

Inserting Eq. (30) into Eq. (26) we find that construc¬ 
tion of OPFs is equivalent to finding a rectangular matrix 
W e U(M, TV) that 


minimizes 


it 


W t t/j l k)t m (k ' b) [/^ k+b) IE 


(31) 
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Here, U\mUA is identified as the enlarged (XI xM) over¬ 
lap matrices projected into the space of M orbitals hj. 

In most cases, the W that minimizes Eq. (31) also sat¬ 
isfies Eq. (28), which then justifies the Taylor expansion 
of uaw ■ However, occasionally this is not the case (for 
example, in strongly covalent materials with a lot of sym¬ 
metry). Therefore, we will introduce a Lagrange multi¬ 
plier A to Eq. (31), which imposes condition Eq. (28). 
With this modification, we now seek matrix W and A at 
a saddle point of the Lagrangian, 


C(W, A) =n 


W'U^m^U { * +b) W 


+Xw ££|[ WU (k)t A (k) lT 
k i— 1 


-1 


(32) 


fin + Hod will be the R = 0 term. Since the R = 0 
term appears only in Hod, it will dominate over Hd for 
an intracell gauge transformation. 

Let us return now back to the optimization problem 
Eq. (26) in question. By construction, the OPFs g 3 ap¬ 
proximately overspan the space of MLWFs near the ori¬ 
gin; in other words, they are related by an intracell gauge 
transformation. Therefore, we are justified in ignoring 
the diagonal part of the spread Hd in Eq. (26). 

With this simplification, the problem of finding ML¬ 
WFs is reduced to finding a rectangular semiunitary ma¬ 
trix W and a real number A which are at a saddle point 
of the Lagrangian, 


C{W, A) =Hi,od 


W^U { ^m^uf +b) W 


For convenience we rescaled the Lagrange multiplier A 
so that A = 1 corresponds to a situation where the rel¬ 
ative importance of the first and second term in the La¬ 
grangian C are equal (w is defined as w = % and 

Wb are k-point weights appearing in the definition of H; 
see Appendix A). 


N 

+A w ££|[ W^A^A^W 

k 2—1 


-1 


(36) 


Inserting here an explicit definition of Hi qd (see Ap¬ 
pendix A) and ignoring the constant term and the l/Aq. 
prefactor, we obtain 


2. Replacing fi with Qi,od 

Now we show that within our approach one can replace, 
in Eq. (26), the total spread H with fli,oD(= Hi + Hod) 
thus ignoring diagonal part of the spread Hd- 

We now examine how the diagonal and off-diagonal 
spread depend on the gauge transformation written in the 
Wannier space. The most general gauge transformation 
of Bloch states is given by Eq. (2) and it involves an 
arbitrary unitary transformation of the states at each k 
point in the Brillouin zone. In the Wannier space, this 
same gauge transformation corresponds to the unitary 
mixtures of WF’s among all unit cells, 

1 0n ) “^ £ £ u ™n I p m) • (33) 

P m 

(P) 

Here the matrix Umn is the Fourier transform of the ma¬ 
trix Umn in Eq. (2). A gauge transformation for which 
(P) 

Umn is nonzero only for P = 0 we will call an intracell 
gauge transformation, since it involves only mixtures of 
the WFs in the same unit cell. 

Let us now start from a set of MLWFs in the home 
cell | On) and see what is the effect of the intracell gauge 
transformation on Hd and Hod- First we will express the 
diagonal and off-diagonal spread in terms of the WFs,' 

Hd = £ I ( R «l r l 0rl )| 2 : (3 4 ) 

n R^O 

Hod = ££ | (Rm|r|On)| 2 . (35) 

m^n R. 

Since the MLWFs are exponentially localized, we expect 
that the dominant term of a gauge dependent spread 


C(W, A) = - % J2 | [w f U^m^ h) U^ +b) W 

k,b 2—1 
N 

+A w ££|[ W t A (k)t A (k) W 


- 1 


k 2=1 


(37) 


Let us now define the following two quantities that are 
independent of W and A: 

M (k - b) = ^ k)t m (k ’ b) 17^ k+b) , (38) 

S (k) = ^(kh^lk) _ Im (39) 

With this simplification, the Lagrangian Eq. (37) now 
simply reads 


c(w,\) = xy a) £ 


(40) 


Here AA“) stands for a collection of M* k,b ) and ma¬ 
trices. The t are the weights associated with the matri¬ 
ces X( a \ with a weight — % for the M( k b ) matrices and 
a weight A w for the S ( k) matrices. Therefore, we have 
reduced a problem of finding MLWFs to the problem of 
codiagonalizing a set of large (M x M) square matrices 
X^ with a single (i.e. k-point independent) rectangular 
(XI xN) matrix W. A mathematically similar approach 
for a square matrix W has been used in Ref. 13 to find 
MLWFs of a localized system. 

In Appendix B we present a numerically efficient algo¬ 
rithm for minimizing Eq. (40), largely following Refs. 14 
and 15. In the following section, we illustrate the OPF 
procedure and empirically validate the approximations 
discussed above. 
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IV. ILLUSTRATION OF OUR APPROACH 

We now illustrate the OPF procedure on a variety of 
systems with chemical bonding ranging from ionic to 
covalent. For predominantly ionic materials we choose 
NaCl, C 1 - 2 O 3 , and LaMnCU. The last two cases have ad¬ 
ditional complexity because they have magnetic and or¬ 
bital order on the transition metals. For predominantly 
covalently bonded materials we choose cubic silicon (c- 
Si), strongly distorted silicon with 20 atoms in the prim¬ 
itive unit cell (Si-20 from Ref. 16), cubic GaAs, and SiC >2 
in the ideal /3-cristobalite structure. 

We computed Bloch wave functions for all seven com¬ 
pounds within the density-functional theory and plane 
wave pseudopotential approach as implemented in the 
Quantum ESPRESSO package. 1 ' The atomic poten¬ 
tials were replaced with ultrasoft 18 pseudopotentials from 
the GBRV 19 library. For the plane wave cutoff, we used 
40 and 200 Ry for the wave functions and charge density, 
respectively. All calculations are done with experimental 
lattice parameters. In the case of Cr 203 we sampled the 
Brillouin zone on a uniform 6 x 6 x 6 k-point grid and 
for all other cases we used a 4 x 4 x 4 grid. 

Using the Bloch wave functions, we computed the over¬ 
lap matrices m( k,b ) between the neighboring Bloch states 
and the overlaps A) k ) between the Bloch states and atom- 
iclike functions that approximately overspan the space of 
MLWFs. For predominantly ionic materials in our test 
(NaCl, Cr 2 0 3 , and LaMn0 3 ), A (k) includes projections 
of Bloch states into all valence atomiclike functions for all 
atoms in the primitive unit cell. For covalently bonded 
materials (c-Si, Si-20, GaAs, and Si0 2 ) some Wannier 
function centers lie on the edge of the primitive unit cell 
(see Sec. Ill), so we included in A^ k) projections onto 
atoms near the edge of the cell. Failing to include these 
additional projections in the case of c-Si yields Wannier 
functions at the computational unit cell boundary with 
spreads two times larger than if we include the additional 
projections. 

We also checked the opposite case by overspanning the 
space of MLWFs even further by including orbitals into 
A (k) that are nominally not in valence (for example, d- 
orbitals in the case of cubic silicon). In this case, the 
final spread for the WFs for the occupied valence band 
complex is unaffected and the matrix elements of W cor¬ 
responding to these additional orbitals is small, as ex¬ 
pected. 

Given matrices m( k,b ) and Af k ' and a choice of the pa¬ 
rameter A we now find matrix W (i.e., OPFs) that mini¬ 
mizes Lagrangian from Eq. (40) using the algorithm de¬ 
scribed in Appendix B. Given W, we construct the u^ v 
to rotate Bloch states into a smooth gauge as described 
in Sec. III. The smoothness of this gauge is quantified by 
first computing the spread fl OPF from the rotated overlap 
matrices in Eq. (26) and then comparing it to the spread 
fi GM at the global minimum. (We define fl GM to be a 
spread of the Wannier functions after running both steps 
of the standard procedure for obtaining MLWFs. For 


TABLE I. Total spread H OPF computed within our approach 
and at the global minimum fi GM for all seven materials stud¬ 
ied. We also give diagonal and off-diagonal components of 
spread in each case (Ho and Hod). The spreads fi OPF are 
obtained using the optimal value of A (see Fig. 1). The units 
for the spreads are A 2 . In the case of Cr 2 03 we wannierize 
only the topmost 12 bands below the Fermi level, and in the 
case of LaMnC >3 we wannierize the top 2 spin-up bands. In 
all other cases we wannierize all valence bands. 




fi OPF 



fi GM 


Total 

Components 

Total 

Components 

Hd 

Hod 

Hd 

Hod 

c-Si 

6.51 

0.00 

0.59 

6.48 

0.00 

0.56 

Si-20 

103.91 

0.05 

14.85 

97.59 

0.04 

8.54 

GaAs 

7.25 

0.02 

0.61 

7.22 

0.01 

0.59 

SiOa 

9.39 

0.00 

1.98 

9.18 

0.00 

1.78 

Cr 2 C >3 

36.04 

0.10 

1.17 

35.74 

0.05 

0.91 

LaMnOs 

14.89 

0.15 

0.17 

14.68 

0.00 

0.11 

NaCl 

4.05 

0.00 

0.80 

4.04 

0.00 

0.79 


convenience, in the first step of finding the global mini¬ 
mum, we do not guess the initial projections but instead 
project into the OPFs obtained from our approach.) 

Figure 1 shows, for all seven cases studied, the ratio 
of the spread U OPF and U GM as a function of A on a 
logarithmic scale. In all cases, the spread U OPF is nearly 
insensitive to the value of A over several orders of mag¬ 
nitude. For example, in the case of GaAs or LaMn 03 
spread U OPF is nearly the same for 0.01 < A < 100. In 
the worst case scenario (c-Si), the spread is still nearly the 
same for 0.1 < A < 2. Therefore, even though in prin¬ 
ciple one may need to vary A to find an optimal value 
of spread, in practice, A ~ 1 is usually a good enough 
choice. 

In each of the seven test cases, the spread U OPF is 
only just 1 % larger than at a global minimum (fl GM ). 
In the worst case situation (Si-20), the spread is only 6 % 
larger than at a global minimum. As mentioned earlier in 
Sec. Ill, this spread could be reduced further by starting 
from OPFs as initial projections and running the second 
step of the standard procedure. 

We give numerical values of fl OPF and U GM in Ta¬ 
ble I along with a decomposition of spread into diago¬ 
nal and off-diagonal components. From here we find an 
additional validation of two simplifications discussed in 
Sec. IIIB. First, Table I shows that linearization of uaw 
is justified since the off-diagonal component of the spread 
bl OPF and U GM is nearly the same. Second, replacing Q 
with Ui od (thus, ignoring diagonal spread) is justified 
within our approach since diagonal spread of U OPF and 
n GM are both very small compared to the total spread. 
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FIG. 1. (Color online) Ratio of fl OPF and fl GM as a function 
of Lagrange multiplier A on a logarithmic scale. 


the WFs are shown on the left, looking down along the c axis. 
The large green dots are La, medium purple dots are Mn, and 
small red dots are O. On the right, we show contour plots of 
the Wannier functions in the plane perpendicular to the c 
axis, cutting through the Mn atom. 


A. Insight gained from the matrix W 

To demonstrate the kind of insight that can be gained 
from analyzing the W matrix, we analyze here in more 
detail case of LaMnC >3 and Cr 2 03 . In both cases, s and 
p orbitals on the neighboring oxygen atoms outside the 
computational unit cell are included in A/ k) in order to 
complete the octahedral coordination of the Cr and Mn 
atoms. 

We studied LaMnC >3 in its low temperature (< 135K) 
A-AFM phase characterized by ferromagnetic ordering 
of the Mn spins in-plane and antiferromagnetic order be¬ 
tween planes. 20 In addition to the magnetic order, Mn d 
states are orbitally ordered, oxygen octahedra are tilted 
and Jahn-Teller distorted. In the following, we focus only 
on the two topmost spin-polarized bands below the Fermi 
level. Analyzing the W matrix we see that the Wannier 
functions for the two topmost bands in LaMnC >3 are dom¬ 
inantly composed of rotated d z 2 components on Mn that 
are oriented perpendicular to each other. This can be 
seen also by analyzing the W matrix for these two WFs, 

|1) « 0.5 |Mnl; d z2 ) + 0.6 |Mnl; d xy ) 

12) « 0.6 |Mn2; d z2 ) - 0.5 |Mn2; d xy ). 

Figure 2 shows a plot of these WFs for the top bands with 
isosurfaces in the left panels and contour plots in the right 
panels. The contours are plotted in the plane perpendic¬ 
ular to the c axis, cutting through the Mn atom. 


Furthermore, The W matrix shows hybridization of 
the Mn d states with the oxygen p states, with the corre¬ 
sponding elements of W having a magnitude of approx¬ 
imately 0.2 (three times smaller than for Mn d states). 
The contribution of the p-like lobes (colored red) can be 
seen in the right panels of Fig. 2 as the large lobes near 
the center. 

Now we analyze the case of Cr 2 03 in more detail. 
Cr 2 C >3 is an antiferromagnetic insulator with four Cr 
atoms in the primitive unit cell. Therefore, we expect 
each Cr 3+ ion to nominally have three occupied d or¬ 
bitals of same spin. These three occupied d orbitals on 
four Cr ions form a complex of 3 x 4 = 12 isolated bands 
that make up the topmost valence bands. Again, ana¬ 
lyzing the W matrix we obtained within our approach 
we find that each of the twelve WFs is a particular lin¬ 
ear combination of all five d orbitals, all having the same 
spin component along the z axis. In fact, there is a large 
degeneracy regarding the particular combination of d or¬ 
bitals that make up the WFs. For example, even slight 
change in A from 1 to 2 gives different linear combina¬ 
tions of d orbitals, while the spread remains nearly the 
same (see Fig. 1). This observation is consistent with the 
fact that the choice of MLWFs is not always unique. This 
was first suggested in Ref. 7 for the case of LiCl. There it 
was found that an arbitrary rotation of the sp 3 orbitals 
on chlorine atoms has no effect on the total spread ft. 




















































V. SUMMARY 


components of the overlap matrices, 


We present an automated procedure for constructing 
maximally localized Wannier functions for an isolated 
group of bands. The extension of our method to the case 
of entangled bands will be the subject of future work. 

Instead of having to guess functions (initial projec¬ 
tions) that approximate the MLWFs as in Ref. 7, our 
approach only requires as input a set of functions that 
overspan the space of MLWFs. In practice, this can 
rather easily be achieved by selecting an appropriate set 
of valence atomiclike functions. 
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Appendix A: Spread functional 


Here we express components of the spread fi, corre¬ 
sponding to A composite bands, as a function of the 
overlap matrices ?rb k ' b ) following Ref. 7, 


^ = at£-K£ 


1 


m. 


(k.b) 


k,b 


On = TW y". Wb ^2 ( -Im hi TO : n b) - b • , (Al) 


k.b 


n ° D = at £ Wb £ 


A k 


m. 


(k.b) 


k.b i^j 


The Wb are the weights of the b vectors connecting neigh¬ 
boring k points (see Sec. 2.1 of Ref. 9), while 


u = vr£ WbbIm ln? 


(k.b) 


A, 


(A2) 


k,b 


We note that the diagonal and off-diagonal parts of 
the spread depend only on the diagonal and off-diagonal 
components of the overlap matrices, respectively. Com¬ 
bining the invariant and off-diagonal parts of the spread 
gives an expression that depends only on the diagonal 


^i,od = Hi + Hqd 


N 




k.b 


i =1 


1 - 


ml 


(k.b) 


(A3) 


Appendix B: Codiagonalization algorithm 


In the main text, the construction of localized Wannier 
functions is recast into the following mathematical prob¬ 
lem. Given a set of large MxAI matrices X^ a \ we wish 
to find a single rectangular semiunitary M x N matrix 
W such that the set of small NxN matrices X^W 
minimize the Lagrangian C , defined in Eq. (40), 


N 

£ <(Q) £|[ w^x^w 


(Bl) 


We parametrize the semiunitary matrix W as follows. 
First, we define W to be first N columns of an M x M 
unitary matrix W. Second, we iteratively parametrize 
the enlarged matrix W as a product (post-multiplication) 
of Givens rotations, 21 

L N M 

W = nn n (B2) 

Z = 1 i = lj=i+l 


Here integer l denotes a particular iteration in the ex¬ 
pansion. 

A Givens rotation R[i,j, 0, </>] is the most general uni¬ 
tary matrix that acts only on ?'-th and j-th rows and 
columns. Therefore, we parametrize R[i,j, 9 , (j)] with two 
angles 9 and <f> as a matrix equal to the identity matrix 
for all elements except for the ii, ij, ji, and jj elements, 


(Ru Ri;j\ _ ( cos 9 e 1 ^ sin <A 

\Rji Rjj) e~’’^sm9 cos 9 ) 

The only diagonal elements of X^ affected by 
R[i,j, 9, (j)\ are X^ and X. Therefore there is no need 
to include in Eq. (B2) cases when both i and j are larger 
than N, since that operation will have no effect on the 
Lagrangian. In addition, we don’t consider cases when 
j < i since that transformation is captured by j > i. 
With this parametrization an arbitrary unitary matrix 
W can be approximated to an arbitrary precision with 
large enough number of iterations, L. 

Let us now see how does a single Givens rotation af¬ 
fect the Lagrangian. For a Givens rotation R[i,j. 9. <[)]. 
the sum of the weighted square moduli of the diagonal 
elements (ii and jj) of a set of rotated matrices R) X^ R 
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are 


15 


with 


E* (a) 


R3X^R 


ll 




R^X {a) R 


2 


jj 


x T Qx + p T x + c 


x T Qx - p T x + c 


(B4) 

A 2 — I 3 

A-! = -2Q 

(B5) 

A 0 = Q 2 - 



where 


The QEP is linearized by introducing 


x T = (cos 29, sin 29 cos </>, sin 29 sin cf >) (B6) 

is a vector with unit norm by construction. The coeffi¬ 
cients of the quadratic forms above [Eqs. (B4) and (B5)] 
depend only on the ii, ij, ji, and jj components of the 
X matrices 



yielding a generalized eigenvalue problem 


(B12) 


(B13) 


z («) z («)t 


q = J 2 t(a) Re 

a 

P =j2 t(a) Re [(Er } + x jf) * z(a) 

a 


(a) 


a: 


(a) 

jj 


(B7) 


where 


z<“> = 1 


y(a) y(a) 

x ii - x jj 

.*{*!?-4 f 


(B8) 


We now consider two cases. First, if j < N both the ii 
and jj diagonal elements enter the Lagrangian C so we 
need to find x that minimizes the sum of Eqs. (B4) and 
(B5), 


E* (a) 

R'X <a) R 


2 + t<“> 

R)X (a) R 


a. 

L 

ii 



33 


= 2x T Qx + 2c. 


(B9) 


This is a quadratic programming problem with the con¬ 
straint that |x| = 1. Here, the Lagrangian is simply 
minimized for x that is the normalized eigenvector corre¬ 
sponding to the minimal eigenvalue of Q. For numerical 
stability, if the first component of x (i.e., cos 2 9) happens 
to be negative we choose —x instead of x. 

In the second case (j > N), only the ii diagonal compo¬ 
nents enters the Lagrangian C so we need to find x that 
minimizes Eq. (B4), 

x T Qx + p T x + c. (BIO) 

The solution of this problem is discussed in Ref. 15 within 
the context of matrix codiagonalization. However we find 
the general quadratic programming solution from Ref. 22 
more numerically stable. Following Ref. 22, we first find 
the minimal eigenvalue x' m in of the quadratic eigenvalue 
problem (QEP) 


(BH) 


Ax = \Bx, 


with 


A = 



A 0 ' 

0 


B = 


'a 2 

0 



(B14) 


(B15) 


This generalized eigenvalue problem we solve using stan¬ 
dard linear algebra techniques. The solution x that min¬ 
imizes Eq. BIO depends on whether Xmin is in the spec¬ 
trum of Q or not. 

If Xmin is not in the spectrum (i.e. not an eigenvalue) 
of Q then the solution is x = {Q — x m ; n I)~ (—p/2). If 
Xmin is an eigenvalue of Q we first define 

U := (Q - XminQ + (-p/2). (B16) 

Here symbol + denotes a matrix pseudoinverse. A non¬ 
trivial solution to Eq. BIO exists only when the following 
conditions are satisfied: 


(Q - XminQu = -p/2 and |u| < 1. (B17) 

Finally, if |u| = 1, then the solution is x = u. Other¬ 
wise (|u| < 1) the solution is x = u + £. Here £ is an 
eigenvector of Q corresponding to Xmin chosen so that 

| C |2 = 1 - |u|2 - 

Once the x is found for a given (i,j) in either of the 
two approaches we determine the corresponding angles 
(9, (j>) from Eq. (B6) and construct the Givens rotation 
R[i, j, 9, </>]. Next we update at each iteration the matrix 
W according to the postmultiplication parametrization 
from Eq. (B2), 


W WR. (B18) 

This iterative procedure over i, j, and l continues until 
the Lagrangian converges. 


{\ 2 A 2 + X^i + Aq)x — 0, 































10 


* jimustafa@berkeley.edu 

1 N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and 
D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012). 

2 X. Wu, O. Dieguez, K. M. Rabe, and D. Vanderbilt, Phys. 
Rev. Lett. 97, 107602 (2006). 

3 J. R. Yates, X. Wang, D. Vanderbilt, and I. Souza, Phys. 
Rev. B 75, 195121 (2007). 

4 F. Giustino, M. L. Cohen, and S. G. Louie, Phys. Rev. B 
76, 165108 (2007). 

5 J. Noffsinger, F. Giustino, B. D. Malone, C.-H. Park, S. G. 
Louie, and M. L. Cohen, Comput. Phys. Commun. 181, 
2140 (2010). 

6 G. Pizzi, D. Volja, B. Kozinsky, M. Fornari, and 
N. Marzari, Comput. Phys. Commun. 185, 422 (2014). 

' N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 
(1997). 

8 I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 
035109 (2001). 

9 A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vander¬ 
bilt, and N. Marzari, Comput. Phys. Commun. 178, 685 
(2008). 

10 A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, 
S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and 
K. A. Persson, APL Mater. 1, 011002 (2013). 

11 K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen, Phys. 
Rev. Lett. 94, 026405 (2005); Phys. Rev. B 72, 125119 
(2005). 

12 R. Sakuma, Phys. Rev. B 87, 235109 (2013). 


13 F. Gygi, J.-L. Fattebert, and E. Schwegler, Comput. Phys. 
Commun. 155, 1 (2003). 

14 J.-F. Cardoso and A. Souloumiac, SIAM J. Matrix Anal, 
and Appl. 17, 161 (1996). 

15 M. Sprensen, L. D. Lathauwer, S. Icart, and L. Deneire, 
Signal Process. 92, 617 (2012). 

16 H. J. Xiang, B. Huang, E. Kan, S.-H. Wei, and X. G. 
Gong, Phys. Rev. Lett. 110, 118702 (2013). 

17 P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, 
C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ- 
cioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fab- 
ris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougous- 
sis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, 
F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, 
L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. 
Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcov- 
itch, J. Phys.: Condens. Matter 21, 395502 (2009). 

18 D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). 

19 K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Van¬ 
derbilt, Comput. Mater. Sci. 81, 446 (2014). 

20 J. B. Elemans, B. V. Laar, K. V. D. Veen, and B. Loopstra, 
J. Solid State Chem. 3, 238 (1971). 

21 G. H. Golub and C. F. Van Loan, Matrix Computations, 
3rd ed. (Johns Hopkins University Press, Baltimore, MD, 
1996). 

22 W. Gander, G. H. Golub, and U. von Matt, Linear Alg. 
Appl. 114-115, 815 (1989). 



